Nonlinear supratransmission in multicomponent systems 
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A method is proposed to solve the challenging problem of determining the supratransmission 
threshold (onset of instability of harmonic boundary driving inside a band gap) in multicomponent 
nonintegrable nonlinear systems. It is successfully applied to the degenerate three- wave resonant 
interaction in a birefringent quadratic medium where the process generates spatial gap solitons. 
No analytic expression is known for this model showing the broad applicability of the method to 
nonlinear systems. 
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Introduction. Nonlinear supratransmission (NST) in 
a medium possessing a natural forbidden band gap is a 
process by which nonlinear structures, gap solitons, are 
generated by an applied periodic boundary condition at 
a frequency in the band gap. Discovered in the pendula 
chain (sine-Gordon model) [1], and further studied for 
fully discrete chain in Refs. [2[ and [ij, it has been ap- 
plied, among others, in Bragg media (coupled mode equa- 
tions in Kerr regime) [4] allowing to explain the experi- 
ments of Ref. [5[ , and also to coupled- wave-guide arrays 
(nonlinear Schrodinger model) pQ- Nonlinear supra- 
transmission results from an instability of the evanescent 
wave profile created by the driving [8[[9| that manifests 
itself above a threshold amplitude. Today this threshold 
has been obtained in single component systems by mak- 
ing use of the explicit solution of the model equation and 
seeking its maximum allowed amplitude at the boundary. 

Predicting the threshold value is of fundamental im- 
portance for physical applications such as soliton gener- 
ation or conception of ultrasensitive detectors. Indeed, 
on the one side NST is a very efficient means to generate 
gap solitons: while an incident single pulse with carrier 
wave at forbidden frequency would be mainly reflected, 
an incident continuous wave excitation easily produces 
gap soliton as experimentally shown in Bragg media [5j. 
On the other side such systems seeded by a CW excita- 
tion slightly below the threshold will be extremely sensi- 
tive to any applied signal, detected either through gen- 
eration of gap solitons or by bistable behavior, see e.g. 

We address in this Letter the practical question of 
evaluating NST thresholds in multicomponent systems 
(where the instability of either wave separately can 
induce soliton formation in all channels), and more- 
over when the system has no explicit solution allow- 
ing for threshold prediction. This is the case with 
second harmonic generation in a birefringent medium 
with quadratic nonlinearity The model is a two- 

component system which does not possess analytic ex- 
pression of solitonlike solutions and which is a key model 
in order to study their existence and stability; see, e.g., 

0. 
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We shall develop a method based on an asymptotic 
solution obtained by asymptotic series expansion, which 
provides an accurate NST threshold prediction. As NST 
requires driving in the forbidden band, the linear evanes- 
cent wave is the natural keystone upon which to build the 
series. The method is restricted to neither the specific 
case of second harmonic generation nor to the quadratic 
nature of the nonlinearity. Moreover, it can be applied to 
a wide class of nonintegrable multicomponent nonlinear 
systems since it does not require known analytical ex- 
pressions for their solutions. The method thus furnishes 
a practical tool highly interesting for further applications 
in any multi-component coupled-wave system. 
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FIG. 1. (color online) Left: intensity plots of a typical numer- 
ical simulation of the D3W system (pQ) with a = 1.1 subject 
to boundary conditions J2} with A = 5, B = 1.3 [the cross 
on the (A, B) plot]. The maximum amplitudes of the emit- 
ted soliton are \ip\ m ~ 4.5 and \<j)\m ~ 5.5. Right: threshold 
curve in the {A, B) plane where dots result from numerical 
simulations [l3[] and the curve is given by the solution of ©. 

Birefringent gap solitons. Let us consider a birefrin- 
gent medium in permanent regime, namely assuming per- 
fect frequency matching. In that case, degenerate spatial 
three- wave model (D3W) reduces to [111] [12] 
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where z) [respectively ^(x, z)e i5z ] is the scaled static 
envelope of the signal wave with frequency uo and wave 
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number k [respectively second harmonic at frequency 
uj' = 2uo and wave number k') and 5 is the missmatch 
wavenumber in the propagation direction z defined by 
k' = 2k — 5. Last, x is the transverse direction and 
a = 1 + S/k f by definition. The system (pQ) is subject to 
the boundary condition 

^(O,*) = Ae 2iz , 0(0, z) = Be iz (2) 

on the strip x > 0, z G [0,P], with vanishing conditions 
as x — >• oo and now with normalization k' = 2 (that is 
a = 1 + 5/2) Q. 

For the sake of simplicity, we only consider here the 
boundary conditions ((2]) , although the wave number may 
in general be chosen different from 1. We see on FigJT] 
the generation, through the evanescent coupling, of bire- 
fringent gap solitons (BGS) formation and propagation 
above a threshold curve in the amplitude plane (A, P). It 
is worth pointing out that the present situation is funda- 
mentally different from studies of initial- value problems 
where field values are given at z = 0, as, e.g., in the the- 
oretical prediction [l5| and experimental realization [l6| 
for nondegenerate 3-wave interaction. 

Asymptotic series solution. Given the boundary val- 
ues (j2j), we seek stationary solutions of the form 

ip(x, z) = u(x)e 2tz , <j)(x, z) = \fot v(x)e lz , (3) 

with real- valued functions u and v vanishing at infinity 
(the factor yfot has been included for convenience). The 
system (pQ), now with S = 2 (a — 1), provides for u and v 
the following parameter free equations 

d 2 x u - 4u + 2v 2 = 0, dlv-v + uv = 0, (4) 

Treating nonlinear terms as perturbative, we first solve 
the linearized equation for v and obtain v = f3e~ x . Sub- 
stituting this result in the equation for u, we find that the 
v 2 term is resonant and generates a solution of the form 
u = (/i + f3 2 x/2)e~ 2x where /i and /3 are two arbitrary 
constants. This is the general solution of the quasilinear 
system d 2 u — Au + 2v 2 = 0, d 2 v = v, that vanishes as 
X — > oo. 

The structure of (|4]) now indicates that u and v may 
be expressed as the following asymptotic series 

oo oo 
n=0 n=0 

Qo(x)=P, P (x)=ii + f3 2 x/2. (5) 

By inspection, the polynomials P n (x) and Q n (x) are of 
degree n (except Po of degree 1) and obey a system 
of differential-recurrence equations obtained by replac- 
ing ([5]) in (UJ). Their coefficients are then recursively 
given in terms of the two independent parameters /? and 
fi. To ensure an accurate determination of the thresh- 
old curve, we have evaluated the first 17 terms of the 



series (j5j), which takes up to a minute on a PC com- 
puter (with MAPLE or MATHEMATICA). The first 
ones are e.g. Q\ = — 3/3 3 /64 - /i/3/8 — (3 3 x/16 and 
Pi = (3 2 ! (24/i + 17/3 2 )/576 + x/3 4 /48. The advantage of 
this method is to be applicable to any system driven in 
a forbidden band, which is essential when no soliton so- 
lution is known. Notice that, in the absence of resonant 
terms in the equations, the polynomials involved in the 
asymptotic series would simply be constants as e.g. in 
the Manakov system below. 

NST threshold prediction. Once the series © has 
been determined up to a given truncation order iV, im- 
posing the boundary conditions u(0) = A and v(0) = B 
leads to the two driving amplitudes A(/3,/j,) and P(/3,/i) 
explicitly expressed in terms of j3 and \i. Assuming that 
the supratransmission threshold curve is given by the 
maximum value of one of those, the other one being 
held constant, we can use a Lagrange parameter A and 
write the extremum condition as d^A — Xd^B = and 
dp A — XdpB = 0. This finally leads to the condition of 
vanishing Jacobian 

J(/3, /i) = {dpA^B) - (d„A){d p B) = 0. (6) 

With a high enough truncation order (N = 16 here), the 
threshold curve is best obtained as the zero contour of 
the surface J(/3, ja) plotted parametrically as a function 
of the amplitudes A(/3, n) and P(/3,/i), as presented in 
Figffl 

The condition ([6]) is symmetric with respect to A and 
P, thus there is no need to specify which maximum am- 
plitude is sought. Moreover condition (J6j) does not de- 
pend on the specific choice of parameters, provided they 
are independent. For instance, using the new parameters 
r] and a defined by 

P = exp(^), / Jj =(a-rj) exp(2r?), (7) 

the NST threshold condition becomes J(a, rj) = 0. The 
interesting point here is that for fixed cr, the new parame- 
ter r] operates a shift of the solutions ([5]) u(x) —> u(x — rj) 
and v(x) — >• v(x — n). Such a translation parameter al- 
ways exists in systems that possess a translation invari- 
ance on the whole x-axis. 

Generalizations. Generalization of the procedure to 
M-component systems is straightforward. To this end, 
let us denote the components by (j) m (x,z) = u rn (x)e lVrnZ 
and their amplitudes by u m (0) — A m ({?7 n }), where r] n are 
the M parameters of the solution (e.g. in the asymptotic 
solution). Finding the NST threshold manifold amounts 
to setting to zero the determinant of the Jacobi matrix 
of the amplitudes with respect to the parameters, that is 

det [Jmn] = 0, Jmn = (8) 

To illustrate our result with another interesting prob- 
lem, we may consider a multicomponent system having 
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soliton solutions, the (integrable) Manakov system [ItJ, 
written here for spatial fields as 

id z ip + d 2 x ip + 2(\<p\ 2 + \ip\ 2 )ip = o, 

«^ + ^0 + 2(H 2 + H 2 )0 = O. (9) 

It possesses the two-parameter soliton solution (ljj 

ijj = e lz sin 6 sech [x — rj) , = e %z cos # sech [x — rj). (10) 

Subject then, on the semi-infinite strip x > and z G 
[0,L], to the boundary condition 

^(0, z) = Ae iz , 0(0, *) = £e i2 , (11) 

and vanishing values as x —> oo, the Manakov system 
possesses solution (flQ|) provided the parameters (77, 0) 
are related to the driving amplitudes A and 5 by 

A = sin sech 77, 5 = cos sech 77 . (12) 

The threshold curve in the (A, B) plane is then obtained 
by solving J(t7, 0) = as given by (|6j). The solution is 77 = 
0, that is, in the variables (A, £?), the circle A 2 -\- B 2 = 1. 
This is illustrated on Figj2] where we display a typical 
soliton formation and the NST threshold curve for which 
the points represent results of numerical simulations [Hj . 




FIG. 2. (color online) Left: typical numerical simulation of 
((9]) submitted to boundary conditions (|lip with A = 0.45 and 
B — 0.9. The maximum amplitudes of the emitted soliton 

are \ip \m ^ 

1.2 and \<fi\m 00 2.3. Right: threshold curve in the 
(A, i^-plane where the dots result from numerical simulations 
131 and the curve is the circle A 2 + B 2 = 1. 

Then one can check that the threshold manifold of 
the M-component Manakov system, idz^m + d 2 (j)m + 
2(5^1^ |^n| 2 )^m = 0, subject to the boundary condi- 
tion m (O, z) = A m e %z , is the M-dimensional sphere 
Si A 2 n = 1. 

Parametric instability and supratransmission. As can 
be seen from FigsQ] and [2j as soon as one of the ampli- 
tudes A or B crosses the threshold curve, the evanescent 
wave profile (soliton tail in integrable models) ceases to 
exist and gap solitons are emitted by the driven bound- 
ary. It turns out that, at least in the models we have in- 
vestigated, the NST threshold manifold also corresponds 



to values of the parameters around which the stability 
of the solutions changes. To substantiate this claim, we 
investigate why the conditions for parametric and NST 
instabilities might indeed be the same by performing lin- 
ear stability analysis of the stationary solutions (|3j) that 
we write in the form 

i){x, z) = (u(x) + U(x, z))e 2iz , (13) 
</>(x,z) = V^(v(x) + V(x,z))e iz . (14) 

where U and V are small perturbations that satisfy, ac- 
cording to (j2j), the boundary condition 

E/(0,z) = 0, V(Q,z) = Q, (15) 

and vanishe at x — >> 00 for all z. Defining the real- valued 
perturbation vectors T r = (U ri V r ) T and I\ = (Ui,Vi) T , 
where index r (i) stands for real (imaginary) part, lin- 
earization of system (pQ) yields 

d z T r + V-Ti = 0, d z Ti - V+T r = 0. (16) 

The matrix differential operators V± are given by 

This is conveniently written as an eigenvalue problem 
by differentiating with respect to z. For the real part 
T r (x, z) we may seek a solution T r (x, z) = ^ u (x) cos(uz) 
and obtain 

V-V+Qu, = co 2 ^ (18) 

with boundary conditions (use d z Ti(0,z) = 0) 

$ w (0)=0, (VM(0) = 0. (19) 

As V± and ^ are real valued, uo 2 is also real valued. 
Then the solution is linearly stable when uo 2 > and un- 
stable for uj 2 < 0, so that marginal instability is reached 
at the bifurcation point 00 = providing the parametric 
instability threshold. 

An essential property of the operator T+ , obtained by 
differentiation of (|4]) with respect to the parameters is 

Thus a two-parameter family of solutions of ([T8]) at uo = 

reads 

for arbitrary constants (a, b) G M 2 . Though it is not the 
most general solution of (fT8]L it seems to be the only one 
able to satisfy the boundary conditions (fT9j) . Requiring 
then <& (0) = 0, with u(0) = A and v(0) = B, eventu- 
ally yields J(t?, a) = 0, namely the parametric instability 
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condition (j6j). Thus the NST threshold condition (j6j), 
i.e., the condition for the maximum allowed amplitudes 
A and B at the boundary x = 0, actually coincides with 
the onset of instability of the solution for the correspond- 
ing critical values of the parameters (here n and a). 

To check this statement, we have computed numeri- 
cally the eigenvalue uj 2 of the differential equation ([T8]) 
around the bifurcation point uj = by varying n and a 
around their critical values r] c and a c defined by the solu- 
tion of J(r],cr) = 0. The result is plotted in Figj3] where 
we have used the asymptotic series solution (u, v) at or- 
der N = 16 as previously. As can be seen from the figure 
marginal instability is actually reached at the criticality 
(f] c and a c ) when uo crosses zero. 
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FIG. 3. Plots of the eigenvalue uo 2 in terms of the parameters 
(77, a) as a function of n — n c for a = a c (left) and function of 
a — a c for 77 = rj c (right) again with a = 1.1. 

Comments and conclusion. The method presented 
here can be readily applied to the simple case of the scalar 
nonlinear Schrodinger equation id z ijj -\- d^t/j -\-2\t/j\ 2 ip = 
to get interesting insight about the occurrence of NST. 
It is found that (i) the asymptotic series solution actu- 
ally sums up exactly to furnish the one-soliton solution, 
(ii) the fundamental parameter is the position n of the 
soliton maximum, (iii) the threshold is indeed the max- 
imum amplitude of this static soliton reached at 77 = 0, 
(iv) the variations of the eigenvalue uj 2 around zero is 
given by uj 2 = — 4r/ + 0(77) that, straightforwardly, gives 
the marginal instability threshold 77 = 0. 

In such single component systems as NLS, the instabil- 
ity occurs always at the maximum amplitude of the solu- 
tion [l9|| . On the contrary, the solution of the D3W model 
does not display any maximum nor any other geometric 
evidence that the NST threshold has been reached. 

In conclusion, we have solved the challenging problem 
of determining the threshold for nonlinear supratrans- 
mission in nonintegrable TV-component systems. This is 
obtained in two steps: first by deriving an asymptotic 
solution based on their linear evanescent profile that de- 
pends on N parameters, and second by solving Eq. (j8]). 
In the parameter space the latter condition results in a 



N — 1 dimensional manifold that determines the change 
of stability of the (asymptotic) solution. Expressed in 
terms of the amplitudes, it gives rise to the sought NST 
threshold. The situation is highly simplified in the case 
of an integrable system, or a system that possesses an ex- 
act static solitonlike solution, since one can work directly 
with the solution to obtain the threshold. 

Finally, since no analytical expression is required, the 
method can be successfully applied to a wide class of 
nonintegrable nonlinear multicomponent systems. 

Work done as part of the programme GDR 3073 
PhoNoMi2 (Photonique Nonlineaire et Milieux Micro- 
structures). 
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